COMPUTING THE GROUND STATE SOLUTION OF 
BOSE-EINSTEIN CONDENSATES 
BY A NORMALIZED GRADIENT FLOW 

WEIZHU BAO * AND QIANG DU t 

Abstract. In this paper, we prove the energy diminishing of a normalized gradient flow which 
provides a mathematical justification of the imaginary time method used in physical literatures to 
compute the ground state solution of Bose- Einstein condensates (BEC). We also investigate the 
energy diminishing property for the discretization of the normalized gradient flow. Two numerical 
methods are proposed for such discretizations: one is the backward Euler centered finite difference 
(BEFD), the other one is an explicit time-splitting sine-spectral (TSSP) method. Energy diminishing 
for BEFD and TSSP for linear case, and monotonicity for BEFD for both linear and nonlinear cases 
are proven. Comparison between the two methods and existing methods, e.g. Crank-Nicolson finite 
difference (CNFD) or forward Euler finite difference (FEFD), shows that BEFD and TSSP are 
much better in terms of preserving energy diminishing property of the normalized gradient flow. 
Numerical results in Id, 2d and 3d with magnetic trap confinement potential, as well as a potential 
of a stirrer corresponding to a far-blue detuned Gaussian laser beam are reported to demonstrate the 
effectiveness of BEFD and TSSP methods. Furthermore we observe that the normalized gradient 
flow can also be applied directly to compute the first excited state solution in BEC when the initial 
data is chosen as an odd function. 

Key words. Bose-Einstein condensate (BEC), Nonlinear Schrodinger equation (NLS), Gross- 
Pitaevskii equation (GPE), Ground state, Normalized gradient flow, Monotone scheme, Energy di- 
minishing, Time-splitting spectral method (TSSP). 
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1. Introduction. Since the first experimental realization of Bose-Einstein con- 
densates (BEC) in dilute weakly interacting gases the nonlinear Schrodinger equation 
(NLS), also called Gross-Pitaevskii equation (GPE) INDIES], has been used extensively 
to describe the single particle properties of BECs. The results obtained by solving 
the NLS showed excellent agreement with most of the experiments (for a review see 
[51 1171 HB| ). In fact, up to now there have been very few experiments in ultracold di- 
lute bosonic gases which could not be described properly by using theoretical methods 
based on the NLS jSlEHj- 

There has been a series of recent studies which deal with the numerical solution of 
the time-independent GPE for ground state and the time-dependent GPE for finding 
the dynamics of a BEC. For numerical solutions of time-dependent GPE, Bao et al. 
[HI E| El E3 presented a time-splitting spectral method, Ruprecht et al. |37] used the 
Crank-Nicolson finite difference method to compute the ground state solution and 
dynamics of GPE, Cerimele et. al. Jl] proposed a particle-inspired scheme. For 
ground state solution of GPE, Edwards et al. presented a Runge-Kutta type method 
and used it to solve Id and 3d with spherical symmetry time-independent GPE |21| . 
Adhikari ^ E] used this approach to get the ground state solution of GPE in 2d 
with radial symmetry. Bao et al. |S] proposed a method by directly minimizing the 
energy functional. Other approaches include an explicit imaginary-time algorithm 
used by Chiofalo et al. |I5j[, a direct inversion in the iterated subspace (DIIS) used by 
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Schneider et al. [2E], and a simple analytical type method proposed by Dodd [T"*"|. In 
fact, one of the fundamental problems in numerical simulation of BEC is to compute 
the ground state solution. 

We consider the NLS equation [J3 HO] 



iipt = - 

y;( x , t) 



- AiP + V(x) V' + /3|V'| 2 V', 



t > 0, 



xe!!C 



0. 



xer = 3fi, t > 0; 



(1.1) 

(1.2) 



where fl is a subset of M. d and V(x) is a real- valued potential whose shape is determined 
by the type of system under investigation, and (3 positive/negative corresponds to the 
defocusing/focusing NLS. is known in BEC as the Gross-Pitaevskii equation 

(GPE) |35| where i\) is the macroscopic wave function of the condensate, t is time, x 
is the spatial coordinate and V(x.) is a trapping potential which usually is harmonic 
and can thus be written as V^(x) = \ + • • • + 7 2 £ 2 i) with 71, • • • , jd > 0. Two 

important invariants of are the normalization of the wave function 



N(iP) 



|^(x,t)| 2 dx 



t > 



and the energy 



i|V^x,t)| 2 + y(x)|,A(x,t)| 2 + ^|^(x,t)| 4 



rfx, t > 0. 



To find a stationary solution of we write 



(1.3) 



(1.4) 



(1.5) 



where fi is the chemical potential of the condensate and <f> a real function independent 
of time. Inserting into i|l.l|) gives the following equation for </>(x) 



M ^(x) = -^A0(x) + V(x) 0(x) + /3|0(x)| 2 0(x), 



x e fi, 



(/>(x) = 0, x g r ; 

under the normalization condition 



/ |0(x)| 2 dx=l. 
Jn 



(1.6) 
(1.7) 

(1.8) 



This is a nonlinear eigenvalue problem under a constraint and any eigenvalue p, can 
be computed from its corresponding cigcnfunction <j> by 



u 



i|V0(x)| 2 + nx)|0(x)| 2 +/3|0(x)| 4 



r/x 



^) + y i^x)! 4 



(1.9) 



The non-rotating Bose-Einstein condensate ground state solution <p g (x) is a real non- 
negative function found by minimizing the energy Ep(cj>) under the constraint (|1.8fl 
[3*2*] . In physical literatures [5J 1131 ITS] , this minimizer was obtained by applying an 
imaginary time (i.e. t — > —it) in ijl.lf) and evolving a normalized gradient flow (see 
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details in the next section). In fact, it is easy to show that the minimizer of Ep((f>) 
under the constraint l|1.8fl is an eigenfunction of l]1.6p . 

The aim of this paper is to prove energy diminishing of the normalized gradient 
flow and present two new numerical methods to discretize the normalized gradient 
flow. This gives a mathematical justification of the imaginary time method which 
is widely used in physical literatures to compute the ground state solution of BEC. 
Energy diminishing of the discretization of the normalized gradient flow is also proven. 
Extensive numerical results are reported to demonstrate the effectiveness of our new 
methods. 

The paper is organized as follows. In section [5] we prove energy diminishing of 
the normalized gradient flow and its discretized version. In section |3 we propose two 
numerical discretizations for the normalized gradient flow. In section 0] numerical 
comparison between the two methods and existing methods, as well as applications 
of the two methods for Id, 2d and 3d ground state solution of BEC, are reported. 
Finally in section [5] some conclusions are drawn. Throughout we adopt the standard 
notation for Sobolev spaces. 

Before we end the introduction, let us note that the NLS is also used in nonlinear 
optics, e.g., to describe the propagation of an intense laser beam through a medium 
with a Kerr nonlinearity |221 140j where ip = '(/'(x, t) describes the electrical field am- 
plitude, t is the spatial coordinate in the direction of propagation, x = (xi, • • • , Xd) T 
is the transverse spatial coordinate and V(x) is determined by the index of refraction. 

2. Normalized gradient flow. In this section we prove energy diminishing of 
a normalized gradient flow and its discretized version. 

2.1. Energy diminishing. Consider the gradient flow 

u t = -Au-V(x)u-/3 \u\ 2 u, t>0, xeft, (2.1) 

tt(x,0) = u (x), xefi, (2.2) 
tt(x,t)=0, xeT, t>0; (2.3) 

where ||u || = 1. Here wc adopt the norm by | • || = || • \\l 2 {Q) an d denote || • = 
II ' llL m (n) with m an integer. Let 

»M) «>a (2.4, 

Then, it is easy to establish the following basic facts: 

Theorem 2.1. Suppose V(x) > for a!lxeS]J>0 and \\u \\ = 1, then 

(i) . \\u{-,t)\\ < |K, 0)|| = IKH = 1 for < t < oo. 

(ii) . For any f3 > 0, 

E p (u(-,t)) < Ep(u(;t')), 0<t'<t<oo. (2.5) 
(Hi). For (3 = 0, 



E {u{-,t)) < E (u(;0)) = E (u ), < t < CO. 



(2.6) 
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Proof: (i). From Ij2.1|l and (|2.3|l . integration by parts, we get 



dt 

= -2 



-HM| 2 = -r / u 2 dx= / 2uu t dx = / 2u 



-Au-V(x)u- (3\u\ 2 u 



dx 



-\Vu\ 2 + V(x)u 2 + /3u 4 



dx < 0, 



< t < oo. 



(2.7) 



This implies the result in (i). 

(ii). From 1)2. ljl . (|2.3Jl and (|1.4(l with ip — u, integration by parts, we get 



dt 



Ep{u) 



-VuVw t + u t (V(-x)u + (3\u\ 2 u) 



dx 



= -2 u t 



-Au- V{x)u~ f3\u\ 2 



2 / \u t \ 2 dx < 0, < t < oo. 



(2.8) 



This implies the result in (ii). 

(hi). From (JOJl with tp = u and /3 = 0, iJHTJ, ||OJ|, H3|l, lE3J and l(2~H|) . 

integration by parts and Schwartz inequality, we obtain 



|Vu| 2 y(x)w 2 



2\\uP || 
Vu • Vut V{~x)u Ut 



dx 



2NI 2 HI 2 

[-\Au + V{*)u] u t 



dx — 


(£w) 


/[ 






Jn . 



|Vu| 



F(x)m 2 



2 u 



"II 
1 



dx-l 5 IMI 



±|Vu| 2 + F(x) u 2 



2 u 



Nl 4 

< , < < < oo. 
This implies l|2~B|) . 



I«H : 



u u t dx - ||u|| \\u t \\ 



(2.9) 
□ 



Remark 2.1. The property $2.5]) is often referred as the energy diminishing 
property of the gradient flow. It is interesting to note that 12. 6)) implies that the 
energy diminishing property is preserved even in the normalized gradient flow when 
[3 = 0, that is, for linear evolution equations. 



Remark 2.2. When /3 > 0, the solution of $U\) - $2~ty) may not preserve the 
normalized energy diminishing property 

Ep{u{;t)) < Ep(u(-,t')), < t! < t < oo. 

In fact, we solve ]2.1\) - ]2.2\) in Id with Q, = R and V(x) = x 2 /2 numerically by the 
time- splitting spectral method ( see details in the next section) for the initial condition 
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uq(x) = e~ x . Figure 0O1 shows, for different [3, the energy Ep(u(-,t)) — 

Ep (u(- , t) / \\u(- , t)\\) under mesh size h = 1/32 and time step k = 0.0001. From the 
figure, we can see that Ep(u) diminishing for < t < oo when fj = 0. But when 
(3 > 0, we have Ep(u) diminishing only for < t < to with some finite to < oo. 




Fig. 2.1. Efj(H) as a function of time in R.emark \2.iA for different /3 (labeled as 5). 

2.2. Normalized gradient flow. Consider the following continuous normalized 
gradient flow 



<f>t = \^<t> - V{x)<t> - [3 |0| 2 + 
<j>(x,t) = 0, xef, 
^(x,0) = 0o (x), xefl 



xeft, i>0, 



(2.10) 

(2.11) 
(2.12) 



In fact, the right hand side of (|2.10|) is the same as (|1.6f> if we view fi^if) as a Lagrange 
multiplier for the constraint l|1.8|) . It readily follows that 



\<K;t)\\- m 



-|V0(x,i)| 2 +F(x)|0| 2 (x,t)+/3|0| 4 (x,<) 



(2.13) 



Furthermore for the above normalized gradient flow, as observed in pfll20|. the solution 
of l|2.10l) also satisfies the following theorem: 

Theorem 2.2. Suppose V(x) > for all x € ft, (3 > and ||</> || = 1. Then 
the normalized gradient flow ^2.10\) - ^2.12\) is normalization conservation and energy 
diminishing, i.e. 



U(-,t)\\= / <^(x,<)dxH|0 o || 2 = i, 

Jet 

^(0) = -2||<M-,;)|| 2 <o, t>o, 



t > o, 



(2.14) 
(2.15) 
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which in turn implies 

Ep^i-M)) > Ef}(<j>(;t 2 )), < ii < t 2 < oo. 



Proof: Multiplying both sides of (|2.10() by (f>, integrating over f2, integration by parts 
and notice (|2~T3T) . J2HH, we obtain 



-A0-V(x)0-/3 3 + M0 (i)0 



' dx 



= 0, 



- 1 V0(x, t) | 2 + U(x)0 2 (x, i) + /30 4 (x, t) 
< > 0. 



djc + ^(t)||0(-,i)|| 



This implies the normalization conservation (|2.14|) . 
Next, direct calculation shows 



= 2 



■ V(f> t + v(*)Wt + P<t?<i>t 



rfx 



-iA0 + y(x)0 + /?0 3 



dx 



= 2/ [-<Mx,i) +/i0(*)'/ , (x,t)]0t dx 

= -2\\M;t)\\ 2 + ^( t ) lX' 0(X ' f)|2dX 
= -2||^(-,<)|| 2 , i>0, 



since £t^(i) is always real and 



j t |j0(x,i)| 2 dx = O 



due to the normalization conservation. Thus, we easily get 

> ^(0(-,t 2 )), 0< t! < t 2 < oo 

for the solution of l|2.10[) . 



(2.16) 



(2.17) 



□ 



Remark 2.3. We see from the above theorem that the energy diminishing property 
is preserved in the continuous dynamic system 1^2. 1(A) . 

Using argument similar to that in |33II39| . we may also get as t — ► oo, <f> approaches 
to a steady state solution which is a critical point of the energy. In non-rotating BEC, 
it has a unique real valued nonnegative ground state solution (fi g (x) > for all x e £1 
|321 - W e choose the initial data </>o(x) > for x S f2, e.g. the ground state solution of 
linear Schrodinger equation with a harmonic oscillator potential |SJ|5]. Under this kind 
of initial data, the ground state solution <fi g and its corresponding chemical potential 
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fj, g can be obtained from the steady state solution of the normalized gradient flow 
<2rnj-f23SJ, i.e. 

<j> g {x) = Urn #x,t), xeO, l*g = l*fi(4>g) = Epfa) + % ( 0*(x)dx. (2.18) 

2.3. Normalized gradient flow via splitting. Various algorithms for com- 
puting the steady state solutions of the normalized gradient flows have been stud- 
ied in the literature. For instance, second order in time discretization scheme that 
preserves the norm normalization and energy diminishing properties were presented 
in 1201- Perhaps one of the more popular technique for dealing with the nor- 
malization constraint is through the following construction: choose a time sequence 
= to < t\ < <2 < • • • < t n < ■ ■ ■ with At n — t n+ \ — t n > and k = max„>o At n . 
To adapt an algorithm for the solution of the usual gradient flow to the case of nor- 
malized gradient flow, it is natural to consider the following splitting (or projection) 
scheme which was widely used in physical literatures [201 023 QH| for computing the 
ground state solution of BEC: 

4>t = ~A<j> - V{y)<t> - H 2 </>, xeO, t n <t<t n+ll n>0, (2.19) 

<f>(x,t) = o, xer, (2.20) 

H0(->*„+l)ll 

0(x,Q) = o (x), xett; (2.22) 

where 0(x, t„) = lim t _^ t ± (f>(x,t) and ||</>o|| = 1. In fact, i|2.19[l is the same as the 
original gradient flow l|2.1fl which can thus be solved via traditional techniques. The 
normalization of the gradient flow is simply achieved by a normalization at each 
discrete time step. 

From Theorem 12. II we get immediately 

Theorem 2.3. Suppose V(x) > for all x <E ft and ||</> || = 1. For = 0, the 
normalized gradient flow is energy diminishing under any time step k and initial data 

E (cf>(;t n+1 )) < Eo(<K; t n )) <■■■< E (<f>(-,Q)) = E (<f> ), n = 0, 1, 2, • • • . (2.23) 
In fact, the normalized step (|2.21|) is equivalent to solve the following ODE exactly 

0t(x,t) = M*(*>fc)<K x >*)> xefi, t n <t<t n+1 , n>0, (2.24) 

0(x,t+) = 0(x,t- +1 ) ) xeSJ; (2.25) 

where 

Ht(t,k) = m(t n+1 ,At n ) = In ||</.(-, t; +1 )|| 2 , t n <t<t n+1 . (2.26) 

Thus the normalized gradient flow can be viewed as a first-order splitting method for 
gradient flow with discontinuous coefficients: 

fa = ^A<t>-V(x)<t>- \cb\ 2 q> + ^(t,k)& xefi, t>0, (2.27) 

<£(x,t)=0, xeT, (2.28) 
0(x,O) = o (x), xe!l. (2.29) 
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Let k — > 0, we see that 



lim fx^t, k) = (j,j,(t) 



■" ' ' l|<K',i)ll 2 J ft 



-|V^(x,t)| 2 +^(x)^(x,i)+/30 4 (x,<) 



dx, 
(2.30) 

which implies that the problem of (|2.27|) - (|2.29|) collapses to (|2.1U[) as k — > 0. 

Remark 2.4. yls we noted earlier, the energy diminishing property in general 
does not hold uniformly for all 4>q and all step size k. Thus, we propose to con- 
sider a modified splitting step which simplifies the computation and yet guarantees the 
monotonicity when it is discretized by BEFD further. 

2.4. Semi-implicit time discretization. To further discretize the equation, 
we here consider the following semi-implicit time discretization scheme: 



1 



k 2 

<T +1 (x) = o, xer, 

0"+ 1 (x) = 0"+ 1 (x)/||0"+ 1 || 



A0 n+i -y(x)(/. n+1 -f3 \ct) n \ 



in|2 ln+l 



x e fi, 



x e Q 



(2.31) 

(2.32) 
(2.33) 



Notice that since the equation i|2.31[l becomes linear, the solution at the new time 
step becomes relatively simple. 
By defining 

V n (x) = V(x)+(3\cf> n ( X )\ 2 , xe!l, 

we may rewrite (|2.31|) as 

ln+1 _ ( m 1 



= -A0 n+1 - V n (x)4> n+1 



(2.34) 



In other words, in each discrete time interval, we may view l|2.31(l as a discretization 
of a linear gradient flow with a modified potential V n (x). 
We now first present the following lemma: 

Lemma 2.4. Suppose V n (x) > for all x € and \\4> n \\ = 1. Then, 



[ \j> n+1 \ 2 dx < f 4> n 4> n+1 dx, 
Jn Jn 



~f l+1 \ 4 dx< / \(f> n \ 2 \4> n+1 \ 2 dx. 



(2.35) 
(2.36) 



Proof: Multiplying both sides of i|2.31[) by 4> n+1 , integrating over SI, and applying 
integration by parts, we obtain 



in i n+1 



dx = —k 



i|V^ +1 | 2 + K(x)|^ +1 | 2 



dx < 
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which leads to lj2.35|l . Similarly, 



In+l|2 



i n+1 - -A0" +1 + kV n (x)4> n+1 



c/x 



\4> 



n+l |2 



6 n+1 | 2 - 2 -0" +1 A0" +1 + 2kV n {-K)\4> n+1 \ 2 



\4> 



n+l |2 



^A0" +1 - kV n (x)4> n+1 



rfx 



= / \<f) n+1 \ 2 |0 n+ T + 3fc|V ( /> n+ T + 2fcK(x)|0" +1 |' 
Jn 



\4> 



n+l |2 



> / |0" +1 | 4 rfx 

JO 



A0" +1 - fcK(x) 



rfx 



This implies <|2~35|> . 



(2.37) 
□ 



Given a linear self-adjoint operator A in a Hilbert space H with inner product 
(•,•), and assume that A is positive definite in the sense that for some positive constant 
c, (u, Au) > c(w, u) for any u G H. We now present a simple lemma: 

Lemma 2.5. For any k > 0, and (/ + /c^w = w, we have 



(u, Au) (v, Av) 
(u,u) ~ (v,v) 



(2.38) 



Proof: Since A is self-adjoint and positive definite, by Holder inequality, we have for 
any p, q > 1 with p + q = pq 



which leads to 



and 



(u,Au) < {u,u) 1/p (u,A?u) 1/q . 
(u,Au) < (u,u) 1/2 (u,A 2 u) 1/2 



(u, Au) (u, A 2 u) < (u, u) (w, A 3 u) . 

Direct calculation then gives 

(u,Au) ((I + kA)u, (I + kA)u) 

= (u, Au) (u, u) + 2k (u, Au) 2 + k 2 (u, Au) (u, A 2 u) 
<{u,Au) (u, u) + 2k (u, u) (u,A 2 u)+k 2 (u,u) (u,A z u) 
= (u, u) ((I + kA)u, A(I + kA)u) . 



(2.39) 

□ 
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Let us defined a modified energy E^n as 



E^(u) = 



|Vu| a + y n (x)|u| s 



dx 



|W| 2 + U(x)M 2 +/3|0"| 2 M 5 



we then get from the above lemma that 

Lemma 2.6. Suppose V(x) > for all x G ft, ft > ant! 



ln+1 



)< 



m+l<l 



/>»+l| 



— Ej,n 



in+1 



/>™+i| 



1 = 1. TTien, 
= £ n(<^ +1 ) <E^{ct> n ) . (2.40) 



Using the inequality l|2.36[l . this in turn implies: 

Lemma 2.7. Suppose V(x) > /or a/Z x e Q, and (3 > 7 t/ien, 



where 



E p {u) 



-\Vu\ 2 + v(x)\u\ 2 + p\ u \ 4 



Remark 2.5. For = 0, the energy diminishing property is preserved in the nor- 
malized gradient flow via splitting \2.iyp - \2.22]) and semi-implicit time discretization 
\2.31\) - y2.33j) . For (3 > 0, we could only justify the energy diminishing on a modified 
energy in two adjacent steps. 

2.5. Discretized normalized gradient flow. Consider a discretization for the 
normalized gradient glow H2.31j) - I|2.33j) ( or a m % discretization of l|2.10|l - l|2.12|l ) 



U 



71+1 



jjn 



k 



-AU 



n+l 



u 



n+l 



u 



n+l 



|{/n+l|| 



n = 0,1,2,- • ■ 



(2.41) 



where U n = «, «§,••• , u^/, & > is time step and A is an (M - 1) X (M - 1) 
symmetric positive definite matrix. We adopt the inner product, norm and energy of 
vectors U — (ui,U2, ■ • • ,«ih) t and V = (v±,V2,-" 7 vm-i) T as 



M-i 



(U,V) = U T V = J2 u 3 v ji \\U\\ 2 = U T U = (U,U), E {U) = U T AU = (U,AU), 
i=i 

(2.42) 

respectively. Using the finite dimensional version of the lemmas given in the previous 
subsection, we have 

Theorem 2.8. Suppose \\U \\ = 1 and A is symmetric positive definite. Then 
the discretized normalized gradient flow \2.4ty is energy diminishing, i. e. 

E (U n+1 ) < E (U n ) <>>-<Eo (U°) , n = 0, 1, 2, • • • . (2.43) 



Furthermore if I + kA is an M -matrix \2Jfl , then {I + kA)~ l is a nonnegative matrix 
(i.e. every entry in it is nonnegative). Thus the flow is monotone, i.e. if U° is a 
non-negative vector, then U n is also a non-negative vector for all n > 0. 

Remark 2.6. If a discretization for the normalized gradient flow \2.31\ - \2.3ty 
reads 
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Suppose B is symmetric and positive definite and p(kB) < 1 where p{B) refers to the 
spectral radius of the matrix B. Then is satisfied by choosing 



A = i ((I - kBy 1 -/) = (!- kBy 1 B. 



Remark 2.7. If a discretization for the normalized gradient flow \2,.31^ - \2.'S'S^ 
reads 

JJ n + l 

U n+1 =BU n , U n+1 = —; , n = 0,1,2,- ■■. (2.45) 



Suppose B is symmetric and positive definite and p{B) < 1. Then \2.43[ ) is satisfied 
by choosing 

A= 1 ~(B-i-l). 



Remark 2.8. If a discretization for the normalized gradient flow \2,.3V^ - \H.'3'3^ 
reads 

fjn+l _ Tjn fj n + 1 

= -BU n+1 - CU n , = -= , n = 0,l,2,---. (2.46) 

k ' 

Suppose B and C are symmetric, positive definite and p{kC) < 1. Then \2.4-3\j is 
satisfied by choosing 

a = (i - key 1 (B + C). 



3. Numerical methods and energy diminishing. In this section, we will 
present two numerical methods to discretize the normalized gradient flow l|2.19[) - 
(I2.22|) . For simplicity of notation we shall introduce the methods for the case of 
one spatial dimension (d = 1) with homogeneous periodic boundary conditions. Gen- 
eralizations to d > 1 are straightforward for tensor product grids and the results 
remain valid without modifications. For d — 1, the problem becomes 

- ' \»-V{x)<f>-p \(p\ 2 (p, x€tt = (a,b), t n <t<t n+1 , n>0, (3.1) 



2 



a 



(x,t 



<f>(x,t n+1 ) = 4>(x,t+ +1 ) = V ' , a<x<b, n>0, (3.2) 

(f)(x,0) = <f> {x), a<x<b, (3.3) 

cf)(a, t) = (f)(b, t) = 0, t>0; (3.4) 

with 

||0o|| 2 = / <j> 2 Q (x) dx = 1. 
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3.1. Numerical methods. We choose the spatial mesh size h = Ax > with 
h = (b — a) /M and M an even positive integer, the time step is given by k — At > 
and define grid points and time steps by 

xj := a + j h, t n := n k, j = 0, 1, • • • , M, n = 0, 1, 2, • • • 

Let (f)™ be the numerical approximation of 4>{xj, t n ) and <j) n the solution vector at time 
t = t n = nk with components iffi. 

Backward Euler finite difference (BEFD) We use backward Euler for time 
discretization and second-order centered finite difference for spatial derivatives. The 
detail scheme is: 



1 



' =^2 iti+i - 2 ^ + - V ( X M ~ P i^f^v i = i, • • • , M - 1, 



k 2h 2 

b*0 = 4>*M = °: 



j = 0,---,M, 7i = 0,l, (3.5) 



$ = fafe), j=0,l,---,M; 
where the norm is defined as 

M-l 
3=1 



Time-splitting sine-spectral method (TSSP) From time t = t n to time 
t = t n+ \, the equation 13.1|l is solved in two steps. One solves 



4>t = ^<t>xx, 

for one time step of length k, followed by solving 

(j> t {x, t) = -V{x)4>(x, t) - /3|0| 2 0(x, t), 



t n <t<t 



n+l, 



(3.6) 



(3.7) 



again for the same time step. Equation l|3.6|l is discretized in space by the sine-spectral 
method and integrated in time exactly. For t e [t n ,t n+ i], multiplying the ODE (|3.7|l 
by 4>(x,t), one obtains with p(x,t) = 4> 2 (x,t) 



p t (x, t) = -2V(x)p(x, t) - 2/3p 2 (x, t), t n <t< t n+1 
The solution of the ODE (|3.8|l can be expressed as 

V{x)p(x,t n ) 



(3.8) 



p(x,t) 



{V(x) + (3p{x, t n )) e 2^(x)(t-t„) _ p p ( Xi tn ) 

P(x,t n ) 

t i + 2pp(x,t n )(t-t n y 



V(x) ± 0, 



V(x) = 0. 



(3.9) 



Combining the splitting step via the standard second-order Strang splitting for solving 
the normalized gradient flow l|3.1(l - (|3.4() . in detail, the steps for obtaining from 
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are given by 



V(xj)e- kv ^ 



V(xj) + /3(1 - e- kv ( x i))\<t><?\ 2 3 
1 



V( Xj ) ± 0, 
V( Xj ) = 0, 



1 ** 



6 j 



M-l 

^ e - fc "</ 2 ^ sin^fo - a)), j = 1, 2, • • • , M - 1, 



V(xj)e~ kV ( x ^ 



V( Xj ) + [3(1 - e~ kV (^))\(j)**\ 2 ^ 3 



1 + W**| 2 



V( Xj ) ? 0, 
V( Xj ) = 0, 



6** 



j = 0,---,M, n = 0, 



where t// are the sine-transform coefficients of a real vector {/ = (uq,ui, 
with uq = um = which are defined as 



M-l 



/i; = — — , E/jr = — ^2 u i sin(/ii(ij - a)), Z = 1,2,- •• ,M- 1 
a i=i 



and 



(^,O) = O (^), j = 0,l,2,...,M. 



(3.10) 

,u M ) T 

(3.11) 



Note that the only time discretization error of TSSP is the splitting error, which is 
second order in k. 



For comparison purposes we review a few other numerical methods which are 
currently used for solving the normalized gradient flow. One is the Crank-Nicolson 
finite difference (CNFD) scheme |23j : 

d>*. — 6 n i 

I 1 2 

- ^ w + m - ™- w+w], i = i , • • • , m - 1 , 

^0 = 0M = °> 
/A* 

J +1 = ^%' J' = 0»---,M, n = 0,l,---, (3.12) 

0° = <M^), i = 0,l,-.-,M. 
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Another one is the forward Euler finite difference (FEFD) method |15j : 
<t>* — <b n 1 

iP~ = ^ y>?+i - 2< t>i + - vfrsw - p i j = i, • • • ,m — i, 

4>l = 4>*m = °. 

^ +1= mV ^°'"-' M ' n = 0,l,--- ) (3.13) 

^°=0o(x,), j=0,l,-.-,M; 

3.2. Energy diminishing. First we analyze the energy diminishing of the dif- 
ferent numerical methods for linear case, i.e. fj — in l|3.1[) . Introducing 



'1- "2- ■•• ) 4>M-l) T ■ 



1 



3=1, 



D = ( d iO(M-i)x(M-i) ' with = 2p S _1 li — 'I = — !> = V" ,M-1, 

~ \ otherwise, 

£ = diag(F(a;i), ^(afc), ••• , 7(ijf_i)) , 

= diag <t> 2 2 , • • • , , with $ = (0!, 2 , • • • , 4> M -i) T , 

G - to)( M -l)x(Af-l) ' Wltn 9]l = T7 2^ 



, , . sin— — -sin— — e fc ^™/ 2 , 
M ^ M M 

m— 1 



F = diag(e- fey ^)/ 2 , e- fcy ^>/ 2 , ••• , e-M^-O/a) . 

Then the BEFD discretization (|3.5[1 (called BEFD normalized flow) with f) = can 
be expressed as 

<j>* — <j>" , $* 

_^_^-(^ + £)$*, n = 0,l,.-.. (3.14) 

The TSSP discretization (|3~TU|) (called TSSP normalized flow) with j3 = can be 

expressed as 

§*** = H<$>** = HG$* =HGH$ n , <$> n+1 = p^|j-, n = 0,l,---. (3.15) 

The CNFD discretization (|3~T2|) (called CNFD normalized flow) with f3 = can be 
expressed as 

^-^ = -l(D+E)$*-±(D+E)$ n , $"+ 1 = ^|l, n = 0,l,---. (3.16) 

The FEFD discretization l|3.13|l (called FEFD normalized flow) with j3 = can be 
expressed as 

= -(£> + £)*", $" +1 = — , n = 0,l,---. (3.17) 

It is easy to see that D and G are symmetric positive definite matrices. Fur- 
thermore D is also an M-matrix and p{D) = (l + cos-^) jh 2 < 2/h 2 and p(G) = 
e -k^l/2 < i Applying Theorem EH and Remarks ^IWTMI^ we have 

Theorem 3.1. Suppose V(x) > and /3 = 0. We have that 
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(i) . The BEFD normalized flow IS. .5)) is energy diminishing and monotone for 
any k > 0. 

(ii) . The TSSP normalized flow L c l.l()\) is energy diminishing for any k > 0. 
(Hi). The CNFD normalized flow is energy diminishing and monotone 

provided that 

2 2h 2 
k ~ 2/h 2 + max., V(xj) ~ 2 + h 2 max., V(xj) ' (3.18) 

(iv). The FEFD normalized flow \3.13\l is energy diminishing and monotone 
provided that 

1 h 2 
k < - = (3 19) 

~ 2/h 2 + maxj V(xj) 2 + h 2 max., V(xj) ' 

For nonlinear case, i.e. (3 > 0, we only analyze the energy between two steps of 
BEFD flow (|3.5() . In this case, consider 

<T>n+l _ fan _ g>n+l 

= -(D + E + /3F($ n ))$ n+1 , $ n + 1 = , (3.20) 

Lemma 3.2. Suppose V(x) > 0, j3 > and ||$"|| = 1. Then for the flow fUfy) . 
we have 

Ep < Ep ($ n ) , E$n < E^ ($") (3.21) 

where 

M-l 



E p ($) = ($,(£> + £ + /JF(*))$) = $ T (D + £)$ + ( 3 - 22 ) 

3=1 

M-l 

($) = ($, (£> + £ + /3F($™))$) = $ T (L> + £)$ + p $ (^™) 2 • ( 3 - 23 ) 

3=1 

Proof: Combining l|3.20|l . (|2.41|) and Theorem we have 

($"+!, (£> + £ + /3F($ n ))<l n+1 
(Z) + £ + /3F($"))$ n+1 ) <^ ? r 



<( ^, (D+E ± pFm m = 

- (<!>™,<I> n ) p 

Similar to the proof of (|2.36(l , we have 

M-l 2 M-l 4 

£(^) 2 (^ +1 ) ^E(^ +1 ) ■ (3-25) 

3=1 3=1 

The required result (|3.21|l is a combination of l|3.25|l . and (|3.24() . □ 
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4. Numerical results. In this section we compare the four different numerical 
discretizations for normalized gradient flow and report numerical results of the ground 
state solutions of BEC in Id, 2d and 3d with magnetic trap confinement potential. We 
also compute the ground state solutions with the potential of a stirrer corresponding a 
far-blue detuned Gaussian laser beam and central vortex state by the methods BEFD 
or TSSP. 

Due to the ground state solution <\> g (x) > for x G ft in non-rotating BEC [221 , in 
our computations, the initial condition (|2.22[) is always chosen such that </>o(x) > and 
decays to zero sufficiently fast as |x| — > oo. We choose an appropriately large interval, 
rectangle and box in Id, 2d and 3d, respectively, to avoid that the homogeneous 
periodic boundary condition Ij3.4(l introduce a significant (aliasing) error relative to 
the whole space problem. To quantify the ground state solution g (x), we define the 
radius mean square 



inns = ||Q!0ff||i*(n) = W J a 2 (j) 2 g (x) dx, a = x, y, or z. (4.1) 

4.1. Comparisons of different methods. Example 1 Normalized gradient 
flow in Id, i.e. d = 1 in H2.19fl - H2.22f> . We consider two cases: 

I. Linear case (f3 — 0) with a double-well potential, 

V{x)= l -{l-x 2 )\ = 0, ^ (x) = — L^e^ 2 / 8 , xel. 

II. Nonlinear case (0 > 0) with a harmonic oscillator potential, 

F(z) = y, = 60, Mx) = -^ UI e- x ^ 2 , xeR. 

The case I is solved on ft = [—16, 16] and the case II on £1 = [—8,8] with mesh 
size h = 35. Figure 4.1 shows the evolution of the energy Ep(<f>) for different time 
step k and different numerical methods. 

From Figure 4.1, the following observations can be made: 

(1) . BEFD is an implicit method and energy diminishing is observed for both 
linear and the nonlinear case under any time step k > 0. The error in the ground 
state solution is only due to the second order spatial discretization. 

(2) . TSSP is an explicit method and energy diminishing is observed for linear 
case under any time step k > 0. For nonlinear case, our numerical experiments show 
that k < guarantees energy diminishing. The error in the ground state solution is 
caused by both the spatial discretization which is spectral accuracy and time splitting 
which is second-order accuracy. From accuracy point of view, large values of k should 
be prohibited. 

(3) . CNFD is an implicit method and FEFD is an explicit method. For both 
schemes, energy diminishing is observed only when the time step k satisfies the con- 
dition (|3.18|) and 13.19JI . respectively. 

To summarize briefly, in general, BEFD is much better than CNFD for computing 
the ground state solution because BEFD is monotone for any k > and CNFD is 
not. TSSP is much better than FEFD. In practice, one can use either BEFD or TSSP. 
BEFD allows the use of much bigger time step k which does not depend on > 0, 
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but the scheme has only second order accuracy in space. At each time step, a linear 
system is solved. In the appendix, we give the detailed BEFD discretization in 2d and 
3d when the potential V(x) and the initial data </>o( x ) have symmetry with/ without 
a central vortex state in the condensate. TSSP is explicit, easy to program, less 
demanding on memory and spectrally accurate in space, but it needs a small time 
step k which depends on the accuracy required and the value of > 0, but not on 
the mesh size h. Based on our numerical experiments given in the next subsection, 
both methods work very well for computing the ground state solution of BEC. 

4.2. Applications to ground state solutions. Example 2 Ground state 
solution of BEC in Id with harmonic oscillator potential, i.e. 

V(x) = ^, ^(aO = -y^e-*"/*, 

The normalized gradient flow (|2.19|) - (|2.22|) with d = 1 is solved on £1 = [—16, 16] with 
mesh size h = | and time step k = 0.001 by using TSSP. The steady state solution 
is reached when max|$ n+1 — $ n < e = 10~ 6 . Figure 4.2 shows the ground state 
solution 4> g (x) and energy evolution for different (3. Table |4~2*1 displays the values of 
4> g (0), radius mean square a; rms , energy Ep{_<t>g) and chemical potential ji g . 

Table 4.1 

Maximum value of the wave function (p g (0), root mean square size a; rms , energy Ep((f>g) and 
ground state chemical potential fi g versus the interaction coefficient (3 in Id. 



a 




^rms 




Ms = M^s) 





0.7511 


0.7071 


0.5000 


0.5000 


3.1371 


0.6463 


0.8949 


1.0441 


1.5272 


12.5484 


0.5301 


1.2435 


2.2330 


3.5986 


31.371 


0.4562 


1.6378 


3.9810 


6.5587 


62.742 


0.4067 


2.0423 


6.2570 


10.384 


156.855 


0.3487 


2.7630 


11.464 


19.083 


313.71 


0.3107 


3.4764 


18.171 


30.279 


627.42 


0.2768 


4.3757 


28.825 


48.063 


1254.8 


0.2467 


5.5073 


45.743 


76.312 



The results in Figure 4.2 and Table 14.21 agree very well with the ground state 
solutions of BEC obtained by a direct minimization the energy functional jH] ■ BEFD 
gives the same results with k = 0.1. 

Example 3 Ground state solution of BEC in 2d. Two cases are considered: 

I. With a harmonic oscillator potential [51 131 \211j. i.e. 

V(x, y) = - {"fix 2 + 7a V) . 

II. With a harmonic oscillator potential and a potential of a stirrer corresponding 
a far-blue detuned Gaussian laser beam \271j which is used to generate vortices in BEC 
123, i.e. 



V(x,y) = t ( 7 > 2 + 7 y) +Wo e- s ^-r°f+y 2 \ 
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The initial condition is chosen as 

For case I, we choose j x = 1, 7 y = 4, wo = S = ro = 0, /3 = 200 and solve the 
problem by TSSP on £1 = [—8,8] x [—4,4] with mesh size h x = g, h y = j§ and time 
step k = 0.001. We get the following results from the ground state solution (fi g : 

x rms = 2.2734, y Ims = 0.6074, 2 (O) = 0.0808, Ep(<t> g ) = 11.1563, fi g = 16.3377. 

For case II, we choose j x = 1, 7 y = 1, wo = 4, S = ro = 1, j3 = 200 and solve the 
problem by TSSP on Q = [—8, 8] 2 with mesh size h = g and time step k — 0.001. We 
get the following results from the ground state solution (j) g : 

x ims = 1.6951, y ims = 1.7144, 2 (O) = 0.034, Ep(<j> g ) = 5.8507, fi g = 8.3269. 

In addition, Figure 4.3 shows surface plots of the ground state solution <fi g . BEFD 
gives similar results with k = 0.1. 

Example 4 Ground state solution of BEC in 3d. Two cases are considered: 

I. With a harmonic oscillator potential f%l H31 \211j. i.e. 

v(x, y,z) = - (i 2 x x 2 + 7 y + 7 2 z 2 ) . 

II. With a harmonic oscillator potential and a potential of a stirrer corresponding 
a far-blue detuned Gaussian laser beam [23 which is used to generate vortex in 
BEC ,23, i.e. 

V(x, y,z)= 1 - ( 7 ^ 2 + 7y V + 7 2 z 2 ) + w e- s ^-^ 2 +y 2 \ 
The initial condition is chosen as 

M*,v,*) = (7a ^I )1/4 ^ feg2+7 ^ 2+7 " 2)/2 - 

For case I, we choose 7x = 1, j y = 2, -f z = 4, w = S = r = 0, (3 — 200 and 
solve the problem by TSSP on Q = [—8,8] x [—6,6] x [—4,4] with mesh size h x = |, 
h y = J^, h z = and time step k — 0.001. The ground state solution cf> g gives: 

x ims = 1.67, y rms = 0.87, z rms = 0.49, 2 (O) = 0.052, Epfa) = 8.33, ii g = 11.03. 

For case II, we choose ^ x = 1, 7^ = 1, ^ z = 2, wo = 4, 6 = ro = 1, /3 = 200 
and solve the problem by TSSP on £1 = [—8, 8] 3 with mesh size h = g and time step 
k = 0.001. The ground state solution cf> g gives: 

x rms = 1.37, y ims = 1.43, z rms = 0.70, </> 2 (0) = 0.025, Epfa) = 5.27, fi g = 6.71. 

Furthermore, Figure 4.4 shows surface plots of the ground state solution <fi^(x, 0, z)M 
BEFD gives similar results with k = 0.1. 
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Example 5 2d central vortex states in BEC, i.e. 

V(x, y) = V(r) = U^+r 2 ), O (*, v) = 0o(r) = r m e^ 2 , < r. 

A \ r* J V7rm! 

The normalized gradient flow is solved in polar coordinate with = [0, 8] with mesh 
size h = J| and time step k = 0.1 by using BEFD (see detail in Appendix A3). Figure 
4.5a shows the ground state solution <f> g (r) with (3 = 200 for different index of the 
central vortex to. Table l4~21 displays the values of <p g (0), radius mean square r rms , 
energy Ep(<fi g ) and chemical potential [i g . 



Table 4.2 

Numerical results for 2d central vortex states in BEC. 



Index to 


<M0) 


'rms 




fig = fJ, {<pg) 


1 


0.0000 


2.4086 


5.8014 


8.2967 


2 


0.0000 


2.5258 


6.3797 


8.7413 


3 


0.0000 


2.6605 


7.0782 


9.3160 


4 


0.0000 


2.8015 


7.8485 


9.9772 


5 


0.0000 


2.9438 


8.6660 


10.6994 


6 


0.0000 


3.0848 


9.5164 


11.4664 



4.3. Application to compute the first excited state . Suppose the eigen- 
functions of the nonlinear eigenvalue problem l|1.6fl . (|l-7f) under the constraint (|1.8fl 
are 

±0 fl (x), ±01 (x), ±02 (x), 

whose energies satisfy 

Ef,(<j> B ) < < Epfa) <■■■ . 

Then (f>j is called as the j-th excited state solution. In fact, <f> g and cpj (j = 1, 2, • • • ) 
are critical points of the energy functional Ep (0) under the constraint 1)1. 8JI . In Id, 
when V(x) = ^- is chosen as the harmonic oscillator potential, the first excited state 
solution 0i (x) is a real odd function, and 4>\{x) = j^t/i x e~ x2 / 2 when (3 = |3i] . 
We observe numerically that the normalized gradient flow (|2.19|l - (|2.22|) and its BEFD 
discretization l|3.5|l can also be applied directly to compute the first excited state 
solution, i.e. 4>i(x), provided that the initial data 4>o(x) in (|2.22(l is chosen as an odd 
function. Here we only present a preliminary numerical example in Id. Extensions to 
2d and 3d are straightforward. 

Example 6 First excited state solution of BEC in Id with a harmonic oscillator 
potential, i.e. 

V(x) = £, ^{ X ) = ^j- A xe-* 2 I\ 

The normalized gradient flow (|2.19|l - (|2.22|) with d = 1 is solved on £1 = [—16, 16] with 
mesh size h = ^ and time step k — 0.1 by using BEFD. Figure 4.5b shows the first 
excited state solution 0i(x) for different j3. Table EOI displays the radius mean square 
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Table 4.3 

Numerical results for the first excited state solution in Id in Example 6. 





x rms 


Ep{4> 9 ) 


Ep(4>x) 


E,3(<p g ) 


Ms 


Mi 


Ms 





1.2247 


0.500 


1.500 


3.000 


0.500 


1.500 


3.000 


3.1371 


1.3165 


1.044 


1.941 


1.859 


1.527 


2.357 


1.544 


12.5484 


1.5441 


2.233 


3.037 


1.360 


3.598 


4.344 


1.207 


31.371 


1.8642 


3.981 


4.743 


1.192 


6.558 


7.279 


1.110 


62.742 


2.2259 


6.257 


6.999 


1.119 


10.38 


11.089 


1.068 


156.855 


2.8973 


11.46 


12.191 


1.063 


19.08 


19.784 


1.037 


313.71 


3.5847 


18.17 


18.889 


1.040 


30.28 


30.969 


1.023 


627.42 


4.4657 


28.82 


29.539 


1.025 


48.06 


48.733 


1.014 


1254.8 


5.5870 


45.74 


46.453 


1.016 


76.31 


76.933 


1.008 



i rms = ||a;0i||L 2 (n)! ground state and first excited state energies Ep((f)g) and Ep(<f>x), 
ratio E/3((j)i)/ E/3((f)g), chemical potentials \i g = fipi^g) and [i\ = ^(</>i), ratio \x\j ' \i g . 

From the results in Table l4~3l and Figure [5p, we can see that the BEFD can be 
applied directly to compute the first excited states in BEC. Furthermore, we have 

lim = lim a = L 

+ Ep{<p g ) fl g 

These results are confirmed with the results in [Hj where the ground and first excited 
states are computed by directly minimizing the energy functional through the finite 
element discretization. 

5. Conclusions. Energy diminishing of a normalized gradient flow and its dis- 
cretization are examined, which provides some mathematical justification of the imag- 
inary time integration method used in physical literatures to compute the ground state 
solution of Bose-Einstein condensation (BEC). Backward Euler centered finite differ- 
ence (BEFD) and time-splitting sine-spectral (TSSP) method are proposed to dis- 
cretize the normalized gradient flow. Comparison between the two proposed methods 
and existing methods shows that BEFD and TSSP are much better for the compu- 
tation of the BEC ground state solution. Numerical results in Id, 2d and 3d with 
different types of potentials used in BEC are reported to demonstrate the effective- 
ness of the BEFD and TSSP methods. Furthermore, extension of the normalized 
gradient flow and its BEFD discretization to compute higher excited states with an 
orthonormalization technique is on-going. 



Appendix: BEFD discretization in BEC when V(x) has symmetry 

In this appendix, we present detailed BEFD discretizations for the normalized 
gradient flows in BEC in 2d and 3d when the potential V(x) and the initial data 
0o(x.) have symmetry with/without a central vortex state in the condensate. Choose 
R > 0, a < b and time step k > with \a\, b, R sufficiently large. Denote the 
mesh size h r = (R — 0)/M and h z = (b — a)/N with M and N two positive integers, 
time steps t n = n k, n = Q, 1, • • • , and grid points Tj = j h r , j = 0, 1, • • • ,M and 
r^i = (j - |) hr, j = 0, 1, • • • ,M + 1, zi = a + lh z , I = 0, 1, • • • ,N. 



Computing the ground state of BEC 



21 



Al. 2d with radial symmetry and 3d with spherical symmetry, i.e. V(x) = V(r) and 
O (x) = o (r) with r = |x| and il = R d with d = 2, 3 in l(2~TO|l - lf2~2^|l . In this case, the 
solution 0(x, t) = 0(r, t) and the normalized gradient flow collapses to a Id problem: 



1 d 

2r d-i 



dr 



V(r)0 - /?|0| 2 0, 0<r<oo, t n <t<t n+1 , (E.l) 



r (O,t) = O, lim 0(r,t) = 0, t > 0, 

r— ► oc 



(r, tn+l) 



< r < oo, n > 0, 



ll<P(->*n+l)ll 

0(r, 0) = o (r) > 0, < r < oo; 
where j|0oj| = 1 and the norm || • || is defined as 



(E.2) 

(E.3) 
(E.4) 



U\\ 2 = C d ^(r,t)r d -'dr 



with 



27T, 
47T, 



c d = 

The BEFD discretization of l|E~T |) -l|E~4 )l is: 

Ci"^ 1 



d = 2, 
d = 3. 



2 hi r d .-\ 



J M- 



i =o, 



,n+l _ 




i-i ~ 


110*11 






J 2 





i = 0,---,M, n = 0,l,--- 



(E.5) 



,M, 



where the norm is defined as 



M 2 



3=1 



A2. 3d with cylindrical symmetry, i.e. V(x) = V(r,z) and 0o(x) = 0o(r, z) with 
r = ^/a; 2 + y 2 and = M d with d = 3 in (|2~H?|) - (|2~2^|) . This is the most popular case 
in the setup of current BEC experiments. In this case, the solution 0(x, t) = 0(r, z, t) 
and the normalized gradient flow collapses to a 2d problem with < r < oo and 

— oo < z < oo: 

~ld_ f df\ d 2 f 
r dr V dr I dz 2 



1 



T/(r,z)0-/3|0| 2 0, t„<t<t n+ i, 



(E.6) 



r (O,z,<) = 0, lim 0(r, z, i) = 0, lim 0(r, z,f) = 0, t > 0, (E.7) 

r — ►oo >±oo 



. A 0(r,^t ra+1 ) 
(r,z,t„+i) = — — — — , n > 0, 

(r,z,O) = o (r,z)>O; 



(E.8) 
(E.9) 
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where ||</>o|| = 1 and the norm || • || is defined as 



/>oo poo 

||0|| 2 = 2tt / / <f) 2 (r,z,t)r dzdr. 



JO J -co 

The BEFD discretization of (lK6ll-(lK9ll is: 



1 

+ 2^! 



1 



2 /i 2 . r, i 

' J 2 



, 2 



j = l,.» ,M-1, l = l,2--- ,N-1, 
**i/ = 0|p 0m_i/ = O, Z = 1, 2, • • • , N - 1, 
^o = ilH = i J = 0,1,...,M. 
4-1 ^i-il 



j = Q,...,M, 1 = 0,1,. ,N, n = Q,l,- 
WW 

Mrj-i,zi), j = l,---,M, l = 0,---,N, 



(E.10) 



.,. / 0. 1.-.-..V. 

where the norm is defined as 



M N-l 



j=l 1=1 



In finding a stationary solution of with a central vortex state, one plugs the 
ansatz 



V>(x,t) 



e -i M t e im9 0( r ^ 



d = 2, 
d = 3, 



r = \/ x 2 + y 2 



into <|1.1[1 instead of l)1.5[l . where m > an integer corresponding to the index of the 
vortex. For more details related to central vortex states in BEC, we refer [181 1291 15H 

A3. 2d central vortex states in BEC, i.e. V(x) = V(r) = \ + r 2 ) and 0o( x ) = 



<j>o{r)with(j)o{Q) = 0, r = yGF+y 2 and Q = R 2 in Q^ - g^ . In this case, the 
solution (/>(x, t) = 0(r, i) and the normalized gradient flow collapses to a Id problem: 

<k = ^^(rTr)-V(r)(f>-0\<t>\*(f>, 0<r<oo, t n <t<t n+1 , (E.ll) 



2r dr \ dr 



0(0,*) =0, 



lim </>(r, t) = 0, i > 0, 

r — >oo 



(r, *n+l) 



A 0(M„+l) 



4>{;t-+i)\\ 
(r,0) = 0o (r) > 0, < r < oo, \v. 



< r < oo, n > 0, 
1 



V irml 



-r 2 /2 



(E.12) 
(E.13) 
(E.14) 
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where 0(0) = 0, \\<po\\ = 1 and the norm || • || is defined as 



r 

||</>|| 2 = 2,t / 
Jo 



4> 2 (r, t)r dr. 



The BEFD discretization of (jEJd) - [jET4j) is: 

<\>* - M l 



2 K r 3 



— r,, i + r.:_ 1 



3+ 



-V(rj) <j>* - p (<f>]Y j = 1, • • • ,M - 1, 

00 = 0M = °: 



6 n+l = 

||0*||' 



i = 0,---,M, 71 = 0, 



(E.15) 



- ru\' 3 

where the norm is defined as 



$ = o (ry), j = 0,l,---,M, 

M-l 

||0*|| 2 = 2^7 r £ r, (0*) 2 . 
A^. 3d central vortex states in BEC, i.e. V(x) = V(r, z) = \ + j 2 r 2 + j^z 2 J 



and 



4>o (x) = (j)o(r, z) with 0o(O, z) = for z e R, 7 r > 0, j z > constants, r = ^ x 2 + y 2 
and = R 3 in (|2.19|) - (|2.22|) . In this case, the solution </>(x, i) = 4>(r,z,t) and the 
normalized gradient flow collapses to a 2d problem with < r < oo and — oo < z < oo: 



i a 



dr \ dr 1 dz 2 



d 2 



V(r,z)(f>- f3cj) 3 , t n <t<t n+1 



<t>(0,z,t) = 0, lim <f>(r, z, t) = 0, lim 0(r,z,<) = O, £ > 0, 

r— »oo 2— >±oo 

0(r, z,t n+1 ) = - - — — — , n > 0, 



ll0(-,^ +1 )|| 
(r ) ! O)=0 o (r,^)>O ! 



1/4 (m+l)/2 ^ 
c .g. = lz y , r m e -(7rr 2 + 7 ^ 2 )/2 



7,3/4^1)1/2 

where </>o(0, z) = for z € R, ||0o|| = 1 an d the norm || • |j is defined as 



(E.16) 
(E.17) 

(E.18) 
(E.19) 



H\\ 2 = 27r 



oo />oc 



<j) 2 {r, z.t)r dzdr. 



The BEFD discretization of (|ET16|l - llET9|) is: 



2 ^ r . 



1 



2 



3 = 1,- •• ,Af- 1, / = 1,2- •• ,N- 1, 
^i = & ( =0, ^0,1, ■■■,7V, <r3* o = 0* M = O, j = 1,1,...,M-1. 



xn+l . 



^, j = 0,---,M, / = 0, ,iV, 77 = 0,1, 



(E.20) 
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where the norm is defined as 

M-l JV-1 

U*\\ 2 = 2nh r h z Y / E(^0 2 ^'- 
j=i 1=1 

The linear system at every time step in Al and A3 can be solved by the Thomas 
algorithm and in A2 and A4 can be solved by Gauss-Seidel iterative method. 
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Fig. 4.1. Energy evolution in Example 1. Left column for case I: a), k = 0.2, c). k = 
e). k = 0.0005. Right column for case II: b). k = 0.05, d). k = 0.01 and f). k = 0.0005. 
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Fig. 4.2. Ground state solution <j> g (x) (labeled as u g ) in Example 2. (a). For (3 = 
0, 3.1371, 12.5484, 31.371, 62.742, 156.855, 313.71, 627.42, 1254.8 (in the order of decreasing 
peak), (b). Energy evolution for different (3 (labeled as 8). 




Fig. 4.3. Surface plots of the ground state solutions <f>g(x,y) (labeled as u 2 ) in Example 3, case 
I (a), and case II (b). 
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0.03. 




Fig. 4.4. Surface plots of the ground state solutions <j>g(x,0,z) (labeled as u 2 ) in Example 4- 
(a). For case I. (b). For case II. 
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(a). r (b). 

Fig. 4.5. (a). 2d central vortex states <j>g( r ) * n Example 5. fi = 200. For m = 1, 2, 3, 4, 5, 6 
(in the order of decreasing peak), (b). First excited state solution <pi(x) (an odd function) in 
Example 6. For [3 = 0, 3.1371, 12.5484, 31.371, 62.742, 156.855, 313.71, 627.42, 1254.8 (in the 
order of decreasing peak). 
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